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ABSTRACT 



Context. Massive stars shape their surrounding medium through the force of their stellar winds, which collide with the 
circumstellar medium. Because the characteristics of these stellar winds vary over the course of the evolution of the star, 
the circumstellar matter becomes a reflection of the stellar evolution and can be used to determine the characteristics 
of the progenitor star. In particular, whenever a fast wind phase follows a slow wind phase, the fast wind sweeps up its 
predecessor in a shell, which is observed as a circumstellar nebula. 
Aims. We make 2-D and 3-D numerical simulations of fast stellar winds sweeping up their slow predecessors to investigate 
whether numerical models of these shells have to be 3-D, or whether 2-D models are sufficient to reproduce the shells 
correctly. 

Methods. We use the MPI-AMRVAC code, using hydrodynamics with optically thin radiative losses included, to make 
O I numerical models of circumstellar shells around massive stars in 2-D and 3-D and compare the results. We focus on 

those situations where a fast Wolf-Rayet star wind sweeps up the slower wind emitted by its predecessor, being either 
a red supergiant or a luminous blue variable. 

Results. As the fast Wolf-Rayet wind expands, it creates a dense shell of swept up material that expands outward, driven 
by the high pressure of the shocked Wolf-Rayet wind. These shells are subject to a fair variety of hydrodynamic-radiative 
instabilities. If the Wolf-Rayet wind is expanding into the wind of a luminous blue variable phase, the instabilities will 
tend to form a fairly small-scale, regular filamentary lattice with thin filaments connecting knotty features. If the Wolf- 
. Rayet wind is sweeping up a red supergiant wind, the instabilities will form larger interconnected structures with less 

regularity. The numerical resolution must be high enough to resolve the compressed, swept-up shell and the evolving 
' instabilities, which otherwise may not even form. 

Conclusions. Our results show that 3-D models, when translated to observed morphologies, give realistic results that 
can be compared directly to observations. The 3-D structure of the nebula will help to distinguish different progenitor 

o 



scenarios. 



■ Key words. Hydrodynamics - Instabilities - Methods: numerical - Stars: circumstellar matter - Stars: massive - Stars: 
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1. Introduction CSM becomes a unique tool to determine the evolutionary 

• ^h path of the central star. A typical example of such an in- 

. As a massive star evolves, it loses a large fraction of its teraction occurs when a massive star (> 30M ) makes the 

r J ; mass in the form of stellar wind. The parameters of the transition from giant star to Wolf-Rayet (WR) star. 
d , wind change during the course of stellar evolution (e.g. 
Ide Jager et al.| [l988), leading to a series of interactions in 

the circumstellar medium (CSM). Whenever the wind pa- As a star reaches the end of its giant phase, it makes 

rameters make a transition from a slow to a fast phase, the a rapid transition from the cool to the hot part of the 

fast wind will collide with its slower predecessor, sweeping Hertzsprung-Russe ll diagram to become a WR star (e.g. 

the older material up in a moving shell. These shells, known lLanger et aD l 19941 ). Thi s transition causes a radica l change 

as circumstellar nebulae, can be observed directly. This will in the wind parameters ([Lamers fc Cassine li3 ll999L and ref- 

be either in emission, if the shells are ionized by UV radi- erences therein). During the giant phase, the wind is rela- 

ation from the central star, or in absorption, as Doppler tively slow (10-200 km s _1 ), due to the low escape velocity, 

shifted absorption features in the spectrum of a separate The WR star, which has a much smaller radius, has a high 

source such as the central star. Since the shells are a direct escape velocity, leading to a fast wind (~ 2000 kms -1 ). As 

result of the changes in wind parameters, which in turn are this fast wind encounters its slow predecessor, it creates a 

caused by the evolution of the star, the morphology of the thin, dense shell, which can be observed as a circumstel- 
lar nebula. Such nebulae, known as WR ring nebulae, ha ve 



Send offprint requests to: A. J. van Marie been found around many WR stars ([Miller fe Chulll993l) . 
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This process is similar to formation of a planetary neb- 
ula, which occurs when a low mass star makes the transi- 
tion from Asymptotic Giant Branch (AGB) to post- AGB 
star ([Kwok et al.| [l978). However, due to the much higher 
kinetic energy of the WR wind as compared to the wind 
of a post-AGB star, the shells of WR nebulae tend to 
be driven by an energy conserving, rather than a momen- 
tum conserving inte raction, lead ing to higher expansion ve- 
locities (chapter 12 lKwokll2000f ). Also, WR stars produce 
very large amounts of high e nergy photons, al l owing them 
to fully ionize their n ebulae ([van Marie et alJ I20Q5L l2QQ7t 
iToala fc Arthurl l201lf ) , rather than the partial ionizatio n 
expected from post-AGB stars (| Garcia- Segur a et aLlll999l ). 

These circumstellar shells, which are pushed by the ther- 
mal pressure of the sho cked WR wind, are subject to linear 
thin-shell instabilities (lVishniadll983l ). Since the material 
in the shell is also much denser than the shocked wind that 
pushes it outward, it can also form Ray leigh- Taylor (RT) 
instabilities. As a result, the structure of these shells can 
become very complicated as has indeed been observed in, 
for e xample, RCW 58 RCW 104, NGC 3199 and NGC 6888 
(E.g.lChu et al.lll983l:lGoudis et al.lfl98 8: Smit h et al.l ll988: 
iDvson fc Ghanbarilll989t iGruendl et al.l l2QQQh . Therefore, 
it becomes necessary to simulate the interaction in more 
than one dimension. 

As a massive star reaches the end of its evolution it 
will die, usually in a spectacular fashion, creating a su- 
pernova and/or gamma-ray burst (GRB). These violent 
events produce fast moving shocks that expand in the CSM. 
The presence of a circumstellar shell can be observed in 
several ways while studying a supernova or GRB. Like 
the central star, the supernova or GRB illuminates the 
CSM with high energy radiation. This ionizes the shell, 
which then in turn becomes visible as electrons and ions 
recombine. This effect i s present in the do uble ring neb- 
ula around SN 1987A ([Burrows et al.lll995l ). which may 
well be the key to understanding the exact nature of its 
proge nitor (e.g. iMorris fc Podsiadlow ski 2005; IChita et~ai1 
2008). The presence of a circumstellar shell can also be ob- 
served as a distinct, blue-shifted absorption f eature in the 
spect r um of some super novae like SN 1998S (Bo wen et al.l 
200(| iFassia et al.l l200lL as its velocity is different from 
the stellar winds. For GRBs this is more difficult to ob- 
serve, as the high energy photons from the GRB and 
its afterglo w will ionize the surr ounding gas to a very 
high degree ([Prochaska et al.ll2007l h though the presence of 
dust in the CSM may compensate for this to some extent 
([Robinson et al.ll2010f ). Alternatively, it may be possible to 
find absorption lines at higher energy levels, if the GRB 
itself is red-shifted sufficiently to bring them within the IR- 
optical part of the spectrum ([Prochaska et al. 2008). 

If an expanding supernova comes in direct contact 
with a circumstellar shell, the interaction can completely 
change the circumstellar environment, depending on the 
ratio of mass between the supernova and the circumstel- 
lar nebula as demon s trated in simulations by, for e x ample , 
[ Tenorio-Tagle et al.l (Il990l Il99lh: iRozvczka etHl (Il993h ; 
ICheyalier fc Dwarkadasjjl995l ): IDwarkadasI ([20051 l2007l ) 
and Ivan Veelen et al.l (l2009t) . This was confirmed observa- 
tionally by |McCravJ ((2005) . When the supernova hits the 
shell, the increased density at the shockfront will create ex- 
tra radiation cau sing a rise in the super nova lightcurve as 
demonstrated by van Marie et al.l ([2010). The collision be- 
tween a supernova and a circumstellar shell was predicted 



in the case of S N1987 A by IChevalier fc Liand (Il989h and 



| Luo fc McCravl (1 1 99 lh. Obse rvations in X-r av dPark et al 
2005), optical ([McCravll2005l ) and infra-rad ([Bouchet et al. 
show brightening due to interaction between the su- 



pernova and the ionized region interior to the circumstellar 
shell, which is denser than the wind. As the supernova rem- 
nant expands further we can expect to observe its interac- 
tion with the shell itself. Whether such an interaction can 
be observed for a GRB is more doubtful, though the transi- 
tion from a 1/r 2 (wind-like) to a constant density medium 
can be visible, depending on th e opening angle of the GRB 
iet (e.g. Ivan Eerten et atll2010D . 

In order to infer the stellar evolution from the 
observable characteristics of the circumstellar environ- 
ment, we need to make detailed numerical models of 
the wind interactions. This has been done pre viousl y 



iFrever et al 



van Marie et al.l ([2005D : 
Dwarkadas] (|2007l ) for the interaction 



Jl996a); 



bet ween a WR wind and a red s u pergiant (RSG) w i nd an d 
bv iGarcia-Segura et all (Il996bh: Ivan Marie et al.l ([20071 ): 
IFrever et al.l (|2003l ): IToala fc Arthurl ([201 ll ) for the inter- 
action between a WR wind and a luminous blue vari- 
able (LBV) wind. Many such simulations were limited 
to 2-D m odels, as was t he la r ge grid of models com- 



puted bv lEldridge et all (|2006l ). IChita et al.l (|2008l ) and 
Ivan Marie et al.l (l2008t) added the effect of rotation and re- 
cently, IToala fc Arthurl (|201lh improved on the earlier 2-D 
models by introducing new physics in the form of radia- 
tive transfer and thermal conduction. Alt hough the effect 
of thermal conduction proved insignificant, Toala & Arthur 
(|201lf ) showed that radiation from the star can play an im- 
portant part in shaping the circumstellar nebula by ionizing 
the shell. These yield a lot of useful insight in the evolution 
of the CSM, but, since they are all limited to 2-D, they 
cannot fully describe what is fundamen tally a 3-D process . 
A single 3-D model was shown in Ivan Marie et all (|2011al ). 
Here we present several simulations at varying resolution 
to investigate the difference between the results of 2-D and 
3-D simulations to determine if and when 3-D models are 

necessary. 

The use of the MPI-AMRVAC code (iM^liani et al.l 120071: 
Ivan der Hoist et al.l 120081 : iKeppens et al.ll2012l and refer- 
ences therein), which can solve the hydrodynamics equa- 
tions for relativistic as well as classical hydrodynamics, also 
helps us prepare for future work, in which the models of the 
CSM can be combined with simulations of supernovae and 
GRBs. As input models we use two stellar wind interac- 
tions, both of which involve a fast wind sweeping up its 
slower predecessor. In the first case, a WR wind sweeps up 
its RSG predecessor. In the second, the WR wind sweeps 
up an earlier LBV wind. For our inpu t models we u s e win d 
parameters based on re c ent work by Bouret et al.l (|2005f ): 
IVink fc de Kotol ([2005D : iMokiem et al.l ([20071) . These re- 
flect the latest developments in our understanding of the 
winds of massive stars. Keeping the physical parameters of 
each model constant, we make both 2-D and 3-D simula- 
tions and vary the resolution of each model. By comparing 
the end results we determine the effect of a change in reso- 
lution and whether 3-D models make a useful contribution 
over the existing 2-D simulations. 

In the following sections, we describe the numerical 
method we use to simulate the interactions (Section [2]) and 
present the results (Section [3]). We discuss the instabili- 
ties that appear in our models in Section SJ Using the 3-D 
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depends on local density, temperature and metallicity. The 
exact calculation of the radiative cooling would actually be: 



R [cm] 



Fig. 1. Density and thermal pressure in the CSM of a WR 
100 years after the end of the RSG phase. This is the start- 
ing point for the simulations labelled with A. 




3.5 o 



Fig. 2. Similar to Fig.[TJ but for the transition from LBV to 
WR. This is the starting point for the simulations labelled 
B. 



models, we show how they would appear to an observer in 
Section [5j Finally, in Section [6] we will analyse and discuss 
the results and in Section [71 present our conclusions. 



2. Numerical setup 

2.1. Governing equations 

We simulate hydrodynamical interactions in the CSM by 
solving the conservation equations for mass: 



| + V.(pu)=0, 



(pu) + V-(puu) 



Vp, 



momentum: 

d_ 
dt 

and energy: 

de _ . . _ . . 
— + V-(eu) + V-(pu) = 



(1) 



(2) 



(3) 



with p the density, u the velocity, p the thermal pressure, 
e the total energy density and mh the hydrogen mass. The 
last equation includes the effect of radiative cooling, which 



— = -n e mA(T), 



(4) 



with n e the free electron particle density and rii the ion 
particle density. Using this equation would require detailed 
knowledge of the local composition and ionization state 
of the gas. For a first approximation we have chosen to 
assume that the gas is fully i onized, which seems a ccept- 
able in the light of results by iToala fc Arthur! (|201lf ) , and 
that hydrogen is the only relevant contributor to the par- 
ticle density, which simplifies the cooling equation to the 
form shown in Eq. (|3]). The cooling function A(T) is a tem- 
perature dependent, which has to be taken from a pre- 
calculated table, taking into account the composition of 
the gas as well as the ionization states at different temper- 
atures. We have chosen to use a ta ble for solar metallic- 
ity from iMellema fc Lundqvistl (|2002f ). which is based on a 
numerical calculation of the energy loss of ionized gas as 
a function of temperature. Using a solar-metallicity based 
cooling curve is justified in that most of the cooling takes 
place in the swept-up shell (due to its high density) of RSG 
or LBV wind material, which can be expected to have the 
same metallicity as the progenitor star had at birth. 

We don't take into accou nt the effect of therma l con- 
duction, but the results from Toal a fc Arthurl (|201lf ) show 
that this is not a significant factor. 



2.2. Basic stellar wind parameters 

As test cases we simulate the interaction between a fast 
WR type star wind and the wind remnant from the 
earlier RSG phase, and the interaction between a WR 
wind and an earlier LBV wind. The RSG wind is slow 
and very dense (V = lOkms -1 , M = 10~ 4 M yr _1 ), 
the LBV wind, while faster (V = 200kms _1 ), has a 
similar high massloss rate (M = 5 x 10 -4 M yr _1 ). 
The WR wind is much faster (2 000kms _1 ), and has 
a lower massloss rate (M = 10 -5 M yr _1 ). These 
particular numbers are b as ed on previous simulation s 
([Garcia-Segura et al.|[l996bl lal: Ivan Marie et al.l2005l I2007D , 
but use a lo wer massloss rate fo r the LBV and WR wind s 
as argued bvlBouret et al.l (|2005f ): IVink fc de Koterl (|2005l ): 
iMokiem et al.l ((2007). For the full parameters of our simu- 
lations, see Table [TJ 

Since the WR star would be a very powerful source 
of high energy radiation, we assume that the CSM 
is fully ionized, leading to a minimum temperature of 
10,000 K throughout the entire computational domain. 
This simplification can be jus tified by ana lytical ap- 
proximation, such as do ne by (Mc Kee et al.l 11984 ) and 
iDvson fc Williams! (|l997l p. 70), which shows that a mas- 
sive star can ionize the surrounding medium out to sev- 
eral par sec, by numeric al results from iFrever et al.l ([2003, 
l2QQ6h : Ivan Marie et al.l ([2005L l2007t ). which incorporated 
a simplified form of p hoto-ionization and recent work by 
IToala fc Arthurl (|201ll ), whic h included radiativ e transfer. 
Observations of NGC 6888bv lMoore et al.l (|2000f ) also show 
that the shell is photo-ionized by radiation from the star. 
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Table 1. Simulation parameters. Simulations labelled with A represent the WR-RSG interaction, simulations labelled 
with B represent the WR-LBV interaction. 



Simulation 


l m (D/y r \ 


71 /T 

J-Vl giant 

[ m o/y r \ 


[km / s] 


^giant 

[km 1 s] 


Domain size 

[PC X A J 


Effective grid 


Al 


1 x 10" 5 


1 x 10" 4 


2000 


10 


5 x 45 


800 x 128 


Bl 


i x icr 5 


5 x 1(T 4 


2000 


200 


5 x 45 


800 x 128 


A2 


1 x 1(T 5 


1 x 10~ 4 


2000 


10 


5 x 45 


1600 x 256 


B2 


1 x 10" 5 


5 x 10" 4 


2000 


200 


5 x 45 


1600 x 256 


A3 


i x icr 5 


1 x 10" 4 


2000 


10 


5 x 45 x 45 


800 x 128 x 128 


B3 


1 x 10" 5 


5 x 10" 4 


2000 


200 


5 x 45 x 45 


800 x 128 x 128 


A4 


1 x 10" 5 


1 x 10" 4 


2000 


10 


5 x 45 x 45 


1600 x 256 x 256 


B4 


1 x 10" 5 


5 x 10~ 4 


2000 


200 


5 x 45 x 45 


1600 x 256 x 256 



2.3. Code 

For our sim ulat ions we use the MPI-AMRV AC code 
dMeliani et all 120071 : Ivan der Hoist et al.l 120081 : 
Kep pens et al.l 120121 and references therein). This is 
a state of the art software package that solves the conser- 
vation equations on grids where the local resolution can 
be changed through adaptive mesh refinement (AMR). 
Together with the extensive use of MPI, this allows us 
to efficiently compute large, multi-D simulations. Both 
the MPI and AMR are absolute necessities for running 
these simulation. E.g. the high-resolution 3-D simulations 
typically take approximately 72-hours, running on 128 pro- 
cessors on either the VIC3 supercomputer at the Flemish 
High Performance Computer Centre or the CINECA 
SP6 in Bologna, Italy. Different physics modules let us 
adapt the code to a specific problem which can include 
magnetic fields and relativistic effects. We have recently 
added a separate modu le to include the effect of op tically 
thin radiative cooling ([van Marie fc Keppensl 1201 lh . this 
modu le uses the exact integration method by iTowns end 
(2009), which improves calculation speed and numerical 
stability. 

2.4. Grid 

We start all simulations in 1-D, using a spherically symmet- 
ric grid. This allows the shock to form without numerical 
instabilities. Our initial 1-D grid has a maximum radius of 
5 pc and a basic grid of 400 points. We allow 7 additional 
levels of refinement, which gives us an effective resolution 
of 51200 points. We initialize the simulation by filling the 
entire grid with wind material according to the parameters 
of the slow (RSG/LBV) wind. At the inner radial boundary 
we allow gas to flow in according to the parameters of the 
fast WR wind. The transition from slow wind to fast wind 
is assumed to be instantaneous, rather than gradual. This 
is a reasonable assumption, since the transition happens on 
the typical free-fall timescale of a star, 



tff K £ = vS' (5) 

with v esc , R* the stellar radius, M* the stellar mass, and 
G Newton's gravitational constant. Even for supergiants, 
this works out to about 10 5 — 10 T s, which is much shorter 
than the timescales on which the nebulae develop (typically 
thousands of years). 



lDwarkad"asl (|2005l [2007) found that starting a simula- 
tion in 1-D influences the final result because turbulence 
in the shocked wind region can influence the shape of the 
wind termination shock. This changes the morphology of 
a RSG shell formed against the shock and does not occur 
in a 1-D simulation. This is not a problem in our case, 
because we do not simulate the formation of such a RSG 
shell, which occurs outside the physical domain of our mod- 
els. Also, we only simulate the first 100 years in 1-D before 
we change to 2-D or 3-D. As the 2-D or 3-D simulation 
starts, we seed the domain with random density variations. 
Thereafter, any side-effects of turbulent shocked wind are 
taken into account. 

After about 100 years, once the shock is formed and 
starts to move outward, we map the 1-D result onto a 2- 
D grid. This approach helps us to avoid num erical prob- 
lems (such as those described by iQuirkl Il994j ). Our grids 
are centred on the equator with a latitudinal opening angle 
of 45° (22.5° above and below the equatorial plane). The 
3-D simulations follow the same pattern with a 45° latitu- 
dinal opening angle. For these simulations the maximum 
radius is again 5 pc. The basic radial grid is 400 points, 
with either 1 or 2 additional levels of refinement. The ba- 
sic angular resolution is 64 gridpoints. (For the resulting 
effective grid resolutions see Table [1]) . 

For our 1-D simulations we use the Total 
Variation Diminishing , Lax-Friedrich (TVDLF) scheme 
(|T6th fc Odst rcil 1996), combined with ' minmod' limiting;. 
For the 2-D and 3-D simulations, we replace the 'minmod' 
limiter with th e more accurate 'van Leer' limiter method 
(|van Leerlll974l ) to compensate for the lower resolution. In 
order to facilitate the formation of instabilities, we ran- 
domly seed both the slow and fast wind with a one-percent 
density variation throughout the entire domain. 



3. Results 

3.1. 1-D models 

The 1-D results, which are the starting points for the 2-D 
and 3-D simulations are shown in Figs. [l]|2j They show the 
free-streaming WR wind, which starts from the central star, 
the shocked wind (starting at 1.6 x 10 17 cm in Fig. [T] and 
6.5 x 10 17 cm in Fig. [2]), the swept up shell (3 x 10 17 cm in 
Fig. [T]and 9.5 x 10 1T cm in Fig. [2j), and the remnant of the 
giant wind. The wind termination shock is near-adiabatic 
in both cases, causing a density jump of approximately a 
factor four. The forward shock (at the front of the shell) 
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Fig. 3. The WR-RSG interaction in low resolution, showing density in cgs of the CSM and the velocity field for simulation 
Al after 7 920 (left panel) and 39 200 years (right panel). At first, the shell clearly shows the development of linear Vishniac 
instabilities. As the shell moves outward, RT instabilities appear as well, leading to large scale disruption of spherical 
symmetry. In the right panel, both Vishniac and RT instabilities remain visible. 




Fig. 4. Similar to Fig. [3j but for simulation A2, which has twice the resolution. Initially, the difference with the low 
resolution run is small. The RT instabilities in the right panel are smaller than their counterparts in Fig. [3] and show 
more structure. Parts of the RT fingers are breaking off and will eventually dissolve in the shocked WR wind bubble that 
is driving the shell. 



is clearly radiative, as the density increases by more than 
three orders of magnitude. 

The thermal pressure is very high in the shocked wind 
region, as it balances the ram pressure of the fast WR 
wind, and remains constant over the contact discontinu- 
ity (at about 2.9 x 10 17 cm in Fig. [Hand 9-5 x 10 17 cm in 



Fig. [2j), since the swept up shell is in pressure equilibrium 
with the shocked wind. 

In the unshocked winds (R < 1.6 x 10 1T cm and R > 
3 x 10 17 cm in Fig. [T] and R < 6.4 x 10 17 cm and R > 
9.3 x 10 17 cm in Fig. [2]) the thermal pressure is compara- 
tively low and is set only by the density and the minimum 
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Fig. 5. 3-D simulation the WR-RSG interaction in low resolution (A3). From left to right: density in cgs of the CSM 
after 7920 and 39 200 years. This figure shows slices through the 3-D grid. These results closely resemble Fig. [3j which 
has the same resolution in 2-D. However, the instabilities in the 3-D model are more pronounced than in 2-D and develop 
quicker. 




0.5 1.0 1.5 2.0 2.5 4.0 5.0 6.0 7.0 8.0 9 .0 10 .0 11 .0 12 .0 13 .0 

(x!0*18 cm) (xlO*18 cm) 



Fig. 6. Similar to Fig. [5j but for simulation A4, which has twice the resolution. As for the low resolution model, the 
instabilities resemble those from the 2-D result, but appear to be further developed. In the later stages they show a 
strong tendency to break off from the shell. 



temperature of 10 000 K. In the WR wind, the ram pressure 
(set by the wind velocity) is about four orders of magni- 
tude higher than the thermal pressure, which means that 
the thermal pressure of the fast wind has no influence on 
the motion of the shell. The thermal pressure in the slow 
(RSG/LBV) wind could become significant, as it has to be 
compared with the forward motion of the shell relative to 



the slow wind, which is much slower than the fast wind 
velocity. Should the thermal pressure of the slow wind gas 
approach the ram pressure felt by the shell as it sweeps up 
the slow wind, then the thermal pressure could influence 
the motion of the shell. As will be discussed in Section [6~T1 
the forward velocity of the shell is about 85 kms -1 for the 
WR-RSG interaction and about 300 kms -1 for the WR- 
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LBV interaction. As a result, in the co- moving frame of 
reference of the shell, the ram pressure exerted by the slow 
wind on the shell (which is actually the effect of the shell 
moving into the slow wind) will be 50 and 800 times larger 
respectively, since the ram pressure increases with the ve- 
locity squared and the sound speed in the unshocked wind 
(which sets the thermal pressure) is about 12 kms -1 . This 
means that the thermal pressure in the slow wind does not 
influence the motion of the shell significantly, which makes 
our assumption of complete ionization of hydrogen in the 
entire grid acceptable. 

3.2. The WR-RSG sequence 

In Fig. [3] we show the results from simulation Al: A 2- 
D model of the expansion of a WR wind driven shell into 
the RSG wind with low resolution. At 7920 years (Fig. [3j 
left panel), the shell already shows instabilities. These are 
thin- shell instabilities of the linear- Vishniac type (Vis hniacl 
Il983l ), caused by ram pressure from the unshocked RSG 
wind on the outside of the shell, balanced by the ther- 
mal pressure of the shocked WR wind. The instabilities 
are small compared to the radius of the shell. As the shell 
moves outward into the slow RSG wind, it forms large scale 
RT instabilities as well as Vishniac instabilities, as a result 
of the density difference between the shocked WR wind and 
the shell it drives. The size of these instabilities is on the 
same order of magnitude as the radius of the shell itself 
(~ 0.5 — lpc). Figure. [3] right panel shows the density at 
39 200 years. Vishniac instabilities are still visible, although 
they are much smaller than the RT instabilities. 

The velocity field shows considerable deviation from 
the general radial motion, as the shocked WR wind flows 
around the RT instabilities. The large RT instability at 
the bottom (right panel of Fig. [3]) blocks the radial mo- 
tion of the shocked wind material completely. This creates 
an "empty" region with lower density (and therefore lower 
thermal pressure) between the instability and the shell. The 
shocked wind material is pushed into this region by the 
thermal pressure gradient, resulting in a strong latitudinal 
motion. Because of the instabilities, there is motion in the 
latitudinal direction inside the shell itself as well. 

The high resolution run (simulation A2) is shown in 
Fig. |U Although generally similar to simulation Al, the 
instabilities tend to be smaller (by about 50 percent) and 
show much more structure, which is understandable, since 
the smaller structures could not be resolved by the low res- 
olution grid. The initial thin-shell instabilities have larger 
amplitudes than in the low resolution model, due to a 
shorter growth time (see Section [4]). Also, parts of the 
larger instabilities are breaking off and dissipating into the 
shocked WR wind, which does not occur in the low resolu- 
tion simulation. 

The results of the low resolution 3-D simulation of the 
WR-RSG interaction (A3) are shown in Fig [5] These have 
the same resolution as Al (Fig. [3j but on a 3-D grid. The 
figures show a 2-D slice through the 3-D grid, which can be 
compared directly to the 2-D results. Clearly, the 2-D and 
3-D simulations produce qualitatively similar results, as the 
location of the shell and termination shock are nearly iden- 
tical. In the first snapshot, the Vishniac instabilities in the 
3-D model are larger than in the 2-D simulation and they 
appear to grow more quickly. It also shows the early stages 
of RT instabilities, which are much more developed then in 



the 2-D model. In the later stage (Fig. [5j right panel), the 
large scale RT instabilities dominate the morphology, which 
closely resembles the 2-D result. Their shape seems to be 
somewhat more longitudinal, but the difference is small. 

The results of the high resolution run (A4) are shown in 
Fig. [6J The instabilities seem to grow more quickly in the 
3-D model than in the 2-D simulation (though the differ- 
ence is less pronounced than for the low resolution model), 
with the early result (left panel of Fig. [6]) showing well 
developed Vishniac and RT instabilities. This is a direct 
result of the difference between 2-D and 3-D geometry. In 
a 2-D model, gas in the instability can only expand in the 
latitudinal direction. In the 3-D model it can move in the 
longitudinal direction as well, giving it an extra degree of 
freedom (See Sect. [6]). The result after 39 2000 years (right 
panel of Fig. [6]), resembles the 2-D result at similar resolu- 
tion. However, the individual instabilities are smaller than 
in 2-D and have a strong tendency to detach from the shell. 

3.3. The WR-LBV sequence 

Figure [7] shows the low resolution result for the WR-LBV 
interaction in simulation Bl after 4 010 and 11800 years. 
Unlike the WR-RSG interaction, the instabilities remain 
small (< 0.1 pc), although both linear- Vishniac and RT in- 
stabilities occur. There are several causes for this effect. 
First, despite its higher massloss rate, the density in the 
LBV wind is lower than in the RSG wind, due to its higher 
velocity (p ex v _1 ). As a result, the density difference be- 
tween the WR wind and the LBV wind is smaller than 
between the WR wind and the RSG wind. This also in- 
fluences the size of the RT instabilities, which depends on 
the density contrast. Second, the shell travels much faster 
(~ 330 kms -1 vs. ~ 85 kms -1 ), giving the instabilities 
much less time to grow. The flow is typically radial through- 
out the grid, with only small deviations around the insta- 
bilities. 

The high resolution WR-LBV simulation (simula- 
tion B2, Fig. [8]) shows considerable difference with the low 
resolution run. The low resolution run clearly has trou- 
ble resolving the instabilities, which appear to be identical 
as they follow the shape of the gridcells. The high resolu- 
tion model resolves the instabilities properly so they show 
local variation. In the high resolution model, one of the 
Vishniac instabilities grows to such a size that it breaks up 
the spherical symmetry of the shell and causes a turbulent 
flow pattern. Apart form this large instability, the flow is 
almost completely radial, with only small local variations. 
The RT instabilities remain thin and don't grow into the 
large clumps that occur in the WR-RSG simulations. This 
is due to the much faster forward motion of the shell, which 
does not give the instabilities time to expand and the lower 
density contrast between the shell and the shocked wind 
region. 

The evolution of the large instability deviates consider- 
ably from that of the shell in general. Its radial expansion is 
much faster and the small RT instabilities at the base tend 
to block part of the flow from the shocked wind region. 
As a result, the pressure driving the instability outwards 
becomes anisotropic and it starts to take on so me of the 
chara cteristics of a non-linear Vishniac instability (Vishniac 
Il994t ). The instability never fully transitions to the non- 
linear regime because the thermal pressure in the shocked 
wind region remains larger than the ram-pressure. The na- 



7 



van Marie & Keppens: Circumstellar shells around evolved massive stars 




Fig. 7. The WR-LBV interaction in low resolution. This figure shows the density in cgs. as well as the velocity field for 
simulation Bl after 4 010 and 11800 years. In the left panel, the shell mainly shows small Vishniac instabilities. In the 
right panel, some evidence of RT instabilities is visible, but unlike for the WR-RSG interaction (Figs.[3]|l]), they are small 
compared to the size of the shell. 




Fig. 8. Similar to Fig. [7J but for high resolution (simulation B2). Not only do the RT instabilities develop much further 
than in the low-resolution model, but one of the Vishniac instabilities grows to such a scale, that it deforms the general 
shape of the shell, with the hot gas of the shocked WR wind breaking out of the shell. This does not occur in the low 
resolution model. 



ture of the forward shock of the instability also changes. 
Due to its high radial velocity shock-heating increases and 
the shock becomes partially adiabatic. This leads to an in- 
crease in the thickness of the shell (Fig. [8j right panel), 
which stops the evolution of the Vishniac instability as the 



curvature of the forward shock starts to decrease relative 
to the triangular shape of the contact discontinuity. 

In order to investigate whether the large distortion of 
the shell that occurs in simulation B2 is representative for 
higher resolution solutions, we repeat the simulation with 
one more level of grid refinement, which gives us an effective 
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Fig. 9. Similar to Fig. [8j but with double the resolution. A break-out similar to the one shown in Fig. [8] occurs in this 
simulation, but it does not grow to the same size and disappears over time. 




Fig. 10. Slices through the 3-D simulation of the WR-LBV interaction in low resolution (B3). From left to right: density 
in cgs of the CSM at 4 010, 7920 and 11 800 years. These results closely resemble Fig. [71 which has the same resolution 
in 2-D. 



resolution of 3200 x 512 gridpoints. It is to be noted that 
the random density variations in fact lead to different re- 
alizations (on top of the resolution differences), preventing 
us to reach true numerically converged states. The higher 
resolution result is shown in Fig. [9] Although a large scale 
instability does form in the earlier stages of the evolution 
of the shell, it does not grow to the same scale as in Fig. [51 
Moreover, the large instability spreads out in the latitudi- 
nal direction, rather than retain its triangular shape. The 



small scale instabilities are similar to those that occur in 
simulation B2. We therefore speculate that the large-scale 
deformations from simulation B2 represent rare events. 

Figures [10] and [TT] show slices through the 3-D simu- 
lations of the WR-LBV interaction (simulations B3 and 
B4). Both simulations mainly resemble the low resolution 
2-D model (simulation Bl, Fig. [7]). The low resolution 3-D 
model shows more variation in the shape of the instabil- 
ities than the 2-D model. Comparisons between the high 
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Fig. 11. Similar to Fig. [TUJ but for simulation B4, which shows the same interaction in high resolution. Compared to 
the 2-D model at the same resolution (simulation B2, Fig. [8]), these results stand out due to the absence of any large 
instability. 



resolution models are difficult due to the absence of the 
single large instability that occurred in the high resolution 

2- D model (Fig. [8]). Since it does not occur anywhere in 
the 3-D models (See also Figs. [T6l and [T7]h such large scale 
structures are probably rare, though a "blow-out" has been 
observed in, for example, RCW 58. A comparison between 
simulation between simulation B4 and the extra-high res- 
olution model shown in Fig. [9] shows that the results are 
very much alike, though, once again, the extra resolution 
of the latter helps to resolve the instabilities. As in 2-D, the 

3- D simulation shows more structured instabilities at high 
resolution, which have a tendency to break off from the 
shell. The instabilities seem to have a tendency to cluster 
together in small groups. These groups are the result of the 
filamentary 3-D structure of the instabilities (See Sect. [5j) 
and denote the places where the slice intersects with the 
knots between filaments. 



4. Instabilities 

There are several different mechanisms for instabilities at 
work in our hydrodynamical models of circumstellar shells. 
Initially, in all models , we see linear thin-shell, or linear 
Vishniac, instabilities (lVishniadll983l ). These occur when- 
ever a thin shell is compressed between thermal pressure 
on one side and ram pressure on the other as is the case 
in all our simulations. They are most visible in the earlier 
stages of the shell formation because of their short growth 
time, which is compa rable to the sou nd crossing time of the 
shell (Eqn. 2.23 from lVishniac1ll983[ ). Since the initial shell 
has a cross-section of less than l/10 th of a parsec and the 
sound speed in the shell (at 10000 K) is about 12 kins" 1 , 
this gives us a formation time of a few thousand years, 
short enough to be visible in the first image (left panels in 
Figs. [3J[8] and [MED) 



The non- lin ear variant of the thin-shell instability 
([Vishniad 1 19941 ) could theoretically occur if the shell was 
subjected to ram pressure from both sides. For this to hap- 
pen, the wind termination shock would have to be radia- 
tive, which is possible at high densities. Because the density 
of the wind is high close to the star, all wind termination 
shocks are initially radiative. However, this is only for a 
very short period of time, as the wind termination shock 
has already become adiabatic after 100 years in the 1-D 
models (see Section ETT]) . This does not give the non-linear 
thin-shell instabilities enough time to grow. 

In the later stages of their evolution, the shells are dom- 
inated by RT instabilities. These result from a situation 
where a low-density gas (the shocked WR wind) acceler- 
ates a much denser gas (the RSG or LBV wind). The size 
of these instabilities is much larger than for the thin-shell 
instabilities, and in the case of the WR-RSG interaction 
they grow to such a scale that they can deform the overall 
shape of the shell. 

The third form of instability is the radiative cooling in- 
stability. This is caused by the strong density dependence 
of the radiative cooling (oc p 2 ). As a result, high density 
regions cool much quicker than low density regions leaving 
them with a lower thermal pressure. This causes them to 
be compressed by the surrounding medium, which in turn 
increases their density. Typically, such instabilities occur on 
a small scale, which makes them difficult to resolve, lead- 
i ng to numerical problems; especially in 2-D simulations 
van Marie fc Keppensll2011l) . They cause local density dif- 
ferences in the fast cooling shell, which can serve as a start- 
ing point for the larger thin-shell and RT instabilities. We 
limit the growth of radiative cooling instabilities by main- 
taining a minimum temperature of 10000 K. This stops the 
cooling clumps from being compressed further, once they 
reach this temperature. 
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Fig. 12. This figure shows the power spectrum calculated from a 1-D Fourier transform along the latitudinal axis for 
the final results of simulations Al and A3, being 2-D vs. 3-D (left panel) and their equivalents at higher resolution, 
A2 and A4(right panel). All simulations show a peak at about 0.8 per latitudinal degree, corresponding to the large 
RT instabilities. The higher resolution simulations (A2 and A4) show more high order maxima than the low resolution 
models. Also, the peaks corresponding to instabilities tend to be higher with respect to the zero-th order peak for the 
3-D simulations than for the 2-D models 




frequency [degree 1 ] frequency [degree 1 ] 



Fig. 13. Similar to Fig.[l2j but for the final results of simulations Bl and B3 (left panel) and B2 and B4 (right panel). As 
for the WR-RSG interaction, the higher resolution models tend to show more higher order maxima and the 3-D models 
have higher peaks corresponding to the instabilities. The power spectrum of simulation B3 shows a peak at about 0.05 
per latitudinal angle, which corresponds to the single large instability. 



A fourth instability, which cannot occur in our models 
because we neglect detailed radiative transfer, is the photo- 
ionization instability. This effect is observed when only part 
of the gas is photo-ionized. The ionized gas, which is much 
hotter than the neutral gas will tend to expand, leaving it 
with a lower density, which in turn is easier to photo- ionize. 
The neutral gas, which is compressed, becomes denser, re- 
combining more quickly and becoming impossible for the 
ionizing photons to penetrate. This eff ect can be important 
in th e structure of planetary nebulae ([Garcfa-Segura et~ai1 
1999). This can also become important for WR nebulae. 



Models by iToala fc Arthurl (|201lh . which include photo- 
ionization through radiative transfer show that, though the 
temperature of the gas in the shells is typically at about 
10 000 K, indicating that our approximation of complete 
ionization is reasonable, the denser clumps are only par- 
tially ionized and cool to a lower temperature. 

4.1. Powerspectra of the instabilities in circumstellar shells 

To make a quantitative comparison between the 2D and 3D 
results, we calculate the Fourier transform of the particle 
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density along the latitudinal grid axis for the final results 
of each simulation 2-D simulation (right panels of Figs. [3} 
[H and [THH]) as well as the slices through the final results 
of the 3-D simulations (right panels of Figs. [5j[6] and [TOi 
[TT]) . This allows us to identify the recurring pattern of the 
instabilities as well as their wavelength. 

The resulting power spectra for the WR-RSG interac- 
tion, normalized to the peak height of the zero-th order 
maximum are shown in Fig. [12] (for details of the Fourier 
quantification, see Appendix [Aj). Apart from the zero-th or- 
der frequency peak, the result of simulation Al (left panel 
of Fig. [T2]) shows a peak at about 0.08 per degree latitude. 
This corresponds to the large scale RT instabilities of which 
there are four over the 45 degree angle of the simulation. 
A second peak, though very faint appears at about 0.15 
per latitudinal degree. The 3-D simulation (A3) shows a 
similar pattern. The peaks lie slightly closer to the zero-th 
order maximum, which is to be expected since the simu- 
lation only shows three major RT instabilities. The peaks 
corresponding to these instabilities are much higher com- 
pared to 2-D simulation. This can be explained by the fact 
that the RT instabilities are longer and thinner, so the pat- 
tern repeats itself over a larger part of the radial domain. 
The Fourier spectra for simulations A2 and A4 (right panel 
of Fig. [T2]) shows a pattern similar to the low resolution 
models. In this case, the frequency of the instabilities is 
higher for the 3-D model than for 2-D. The 3-D model now 
shows several higher order peaks, which are much higher 
compared to the zero-th maximum. As for the low reso- 
lution models, the repetitive pattern of the instabilities is 
more pronounced in 3-D than in 2-D. 

The results for the WR-LBV interaction are shown in 
Fig. [T3j Both 2-D simulations show peaks at a frequency 
of about 0.3 per latitudinal angle, corresponding to about 
12 waves over the 45° angle of the domain. This frequency 
corresponds with the instabilities noted in Sect. [3) A second 
peak at about 0.5 per latitudinal angle is faintly visible for 
simulation Bl and very clearly for simulation B2. This cor- 
responds to the smaller subdivisions of the initial instabili- 
ties, which are clearly visible for simulation B2. Simulation 
B2 also shows a peak at 0.05 per latitudinal angle. This is 
equal to about 2 waves over the 45° degree domain, which 
corresponds to the shape of the large Vishniac instability. 
The 3-D simulations show a different pattern. Both have 
maxima at about 0.1-0.15 per latitudinal angle and a series 
of peaks at higher frequencies. The peaks at about 0.1-0.15 
are especially pronounced for the low resolution simulation 
(B3). This corresponds to approximately 5 waves over the 
45° degree domain. A close inspection of the density plot 
(right panel of Fig. [TO]) shows that the small RT fingers 
tend to cluster together where the slices intersect with the 
knots in the filamentary structure of the shell (See Sect. [5]), 
forming 4-5 groups. This pattern is less clear for the higher 
resolution model (B4), although the tendency of the fingers 
to appear close together is still visible. 

A similar analysis could be done on observational results 
to identify the patterns of the instabilities and compare 
them to observational results. 

5. Observable 3-D structure 

Using our 3-D simulations, we attempt to simulate how the 
nebulae we produced would appear to an observer. Because 
recombination is the main source of photon production in 



these emission nebulae and the recombination rate is pro- 
portional to the particle density squared, we calculate the 
square of the i on density (n 2 = p 2 /mj) as a measure of 
emissivity (see iGarcfa-Segura et all 1 19991 ) and use it as a 
basis of a volume rendering of the results shown in the 
right panels of Figs [5j [6j [10] and [Til The result is shown 
in Figs. rmfTTl From these figures it becomes clear that the 
nebula will not appear as a smooth shell. Rather, it presents 
itself as a lattice work of high density filaments with empty 
space in between. Of course, the filaments are connected 
through a shell, as can be seen in the slices through the 3- 
D data stack (Figs.lMTT]). However, this shell is so thin, that 
it does not contribute significantly to the emissivity. Only 
the instabilities, which have a high column depth along the 
line of sight become visible. 

Typical values for in the shell are between 10 2 and 
10 5 . To get actual numbers for the emissivity this would 
have to be multiplied with the local value of A(T) (see 
Section I2TT]) . The temperature in the shell is typically of 
the order 10 4 , but since the shock is strongly radiative, 
the radiative cooling causes a significant drop in the gas 
temperature, causing energy gained through the collision 
to be radiated away on a timescale shorter than the hy- 
drodynamical timescale of the gas. As a result, the tem- 
peratures found in the results, which are equilibrium tem- 
peratures, may underestimate the actual radiative temper- 
ature, which should be close to the forward shock temper- 
ature (10 5 - 10 6 K). Typical cooling curves show that solar 
metallicity (a reasonable approximation for the material 
in the shell, which consists of either LBV or RSG wind), 
would give A (T) ~ 10~ 22 — 10~ 21 erg cm 3 s -1 for these tem - 
perat u res (iDalgarno fc McCravl|l97 2; Mac Donald fc Baile 
1 19811: iMellema fc Lundavistl BOOa ISmith et all l200l» 
ISchure et alJl2QQ9h - For our simulations this would result in 



shell-emissivities of the order of 10 _2U — 10 -16 erg cm -3 s _x . 
These estimates give an approximation for emissivity in 
the continuum. In order to find line emissivity, one has to 
use the strength of the emission line, which depends on 
chemical composition and degree of ionization as well as 
temperature and density. Because the shell consists of RSG 
wind material, which is made up of the stellar envelope, the 
chemical composition can be approximated by the ZAMS 
composition of the progenitor star. (A more accurate esti- 
mate can be made based on a stellar evolution model that 
takes into account chemical mixing due to convection.) The 
degree of ionization is a more complicated issue, as it de- 
pends not only on local gas properties, but on t he radiation 
field of the central star (jToala fc Arthur]l201ll ). 

Figures fT4lfT7l show a large difference between the WR- 
RSG interaction and the WR-LBV interaction. The WR- 
LBV shell appears to have a very regular morphology with 
little difference between the local structures. In Fig. [14] in- 
stabilities are approximately the same size, with filaments 
having a width of about 1/100^ of a parsec and distances 
of about l/ltfh of a parsec between knots. In Fig. [T71 the 
number of filaments is much larger. The size of the indi- 
vidual filaments remains very regular, but the distance be- 
tween knots has become smaller. The WR-RSG shell on the 
other hand, is highly irregular, with large scale instabilities 
as well as small filaments. Here the width of the filaments 
varies between l/100 t/l and more than l/10 t/l of a parsec 
and they can be more than a parsec in length. 

The effect of a higher resolution, already shown in 
Section [3] is again visible, though the effect is quantitative 
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Fig. 14. Particle density squared [1/cm 6 ] volume rendering of the low resolution WR-RSG interaction (simulation A3) 
at the same time as the last panel of Fig. [5j 



rather than qualitative. For the WR-RSG interaction, the 
low resolution model (simulation A3, Fig. [T^|) shows less 
fine structure than the high resolution model (simulation 
A4, Fig. [T5]h but the morphological trend is clear. For the 
WR-LBV interaction, the general, fairly regular, structure 
is the same, but in the high resolution model (simulation 
B4, Fig.[T7|) the individual structures are smaller and a sec- 
ondary network of extremely thin filaments has appeared 
that could not be fully resolved by the low resolution model 
(simulation B3, Fig. [T6|) . This corresponds to the second 
peak noted in the Fourier spectrum of simulation B2. 

When comparing these figures to observational re- 
sults, the most obvious comparison can be made with 
the "Crescent nebula", NGC 6888 o f which many obser- 
vations exist in high resolution (e.g. iGruendl et al.l [2000; 
iMoore et al.l l2QQ(5h . This nebula, which is thought to be 
the result of a WR-RSG interaction shows filaments on 
a variety of scales (ranging from l/100 £/l of a parsec to 
scales close to a parsec). Particle densities of clumps in 
NGC 6888 vary considerably, but Hubble Space Telescope 
(HST) obs ervations show clu mps with densities of about 
10 3 cm -3 (M oore et al.l l2000h , which coincides with our 
models, which show particle densities (p/rrih) between 10 2 
and 10 4 . 

5.1. Observed morphology at a distance. 

When looking at the 3-D images in Figs. [T4l[T71 one should 
keep in mind the scale of the colour table. This varies over 



ten orders of magnitude for each figure and the shells have a 
variation of more than three orders of magnitude. For these 
figures, volume rendering has been used to show the entire 
shell. When observing these shells from a distance, only the 
brightest parts of the shell would be visible, leaving more 
empty space between the filaments. In order to approximate 
this result, we have made a series of images [TBI and [T9l in 
which we show only part of the gas, leaving out those parts 
where the density becomes too low (and therefore too faint). 
For this we only use the high resolution simulations (A4 and 
B4). 

Figure [18] shows the frontal aspect of the same shell as 
in Fig. [15] with a lower limit of log(n 2 ) > 2 (left panel), 
log(n 2 ) > 3 (center panel), and log(n 2 ) > 4 (right panel). 
Where originally filaments took up about half of the surface 
area, the open space now starts to dominate as the larger 
filaments appear thinner and the smaller filaments disap- 
pear completely, reducing the visible high density features 
to about ten percent of the surface area. 

For the WR-LBV interaction (simulation B4), we have 
to shift the cut-off values, since this nebula is fainter to 
begin with due to lower densities. The results is shown in 
Fig. [19] with cut-off values log(n 2 ) > 1, 2, and 3. Starting 
out an order of magnitude fainter than its WR-RSG coun- 
terpart, the qualitative structure of this nebula does not 
change much as it becomes fainter, though as in the WR- 
RSG nebula, the size of the voids between the filaments 
increases as the lower density parts of the filaments fade 
away. These voids are not actually empty, but represent 
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Fig. 15. Particle density squared [1/cm 6 ] volume rendering of the high resolution WR-RSG interaction (simulation A4) 
at the same time as the last panel of Fig. [6j 



those parts of the shell, which have no large RT instability. 
Because of this, the cross-section of the high-density region 
is very short, leading to a low column density. Initially the 
high density filaments covered about ten percent of the sur- 
face area. In the final panel, they are reduced to about two 
percent. 

This analysis does not take into account the limits of 
spatial resolution. For example, the HST can resolve struc- 
tures down to approximately 0.1". This means that the 
large filaments of the WR-RSG interaction can be resolved 
for distances of up to about 2 Mpc. On the other hand, the 
smallest filaments in the WR-LBV interaction can only be 
resolved at distances of less than 0.02 Mpc. The Herschel 
Space Observatory infra-red satellite has a spatial resolu- 
tion that is about a factor 350 less than HST and can only 
resolve the small structures at less than 60 pc. Even the 
larger filaments in the WR-LBV shell would appear would 
be impossible for Herschel to resolve at more than 600 pc. 
Therefore, the shells that result from a WR-LBV wind in- 
teraction would appear smooth under most circumstances. 

5.2. Alternative observations 

With the launch of satellites like the Spitzer Space Telescope 
and the Herschell Space Observatory it has become possible 
to observe the morphology of circumstellar nebulae in the 
infra-red. E.g. Spitzer observations, usi ng the Multiband 
Photo-Imager (MIPS) (iRieke et al.ll2004ft. show circumste l- 
lar shells around Wn8 and Wn9-h ([Mauerhan et al.ll2010l ). 



The main contributor to infra-red emission in circumstel- 
lar nebulae is dust. Therefore, in order to predict how our 
models would look in the infra-red we would have to in- 
clude the presence of dust in our si mulations. Even without 
this we can make som e predictions. Ivan Marie et aT1([2011cI ) 
and I Cox et al.l (|2012f ) showed that small dust grains, which 
are by far the most numerous, are tightly bound to the 
gas through the drag force. Therefore, assuming that the 
dust was initially distributed evenly throughout the RSG 
wind, we can assume that the dust density is directly cor- 
related to the gas density and will show the same morphol- 
ogy. This assumption becomes invalid, if significant dust 
formation takes place in the swept-up shell. Should that 
be the case, the newly formed dust would be concentrated 
in the high density blobs, because dust formation depends 
strongly on gas density. Dust emission scales with the dust 
density rather than the density squared, which will make 
the filamentary structure less prominent. Similarly, absorp- 
tion scales with the density, so any observations that show 
a circumstellar shell in absorption will show only limited 
density contrast. Still, since the 3-D plots show differences 
of 5 orders of magnitude or more for n 2 , any observation 
that scales with n will still show a contrast of several orders 
of magnitude. 

The large contrasts in column density limit the infor- 
mation on a circumstellar shell that can be gained from 
the analysis of the spectrum of the central star, since the 
amount of matter observed can change by several orders of 
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Fig. 16. Particle density squared [1/cm 6 ] volume rendering of the low resolution WR-LBV interaction (simulation B3) 
at the same time as the last panel of Fig. [TUJ 



magnitude depending on whether the line of sight passes 
through one of the filaments. 



5.3. The large scale structure of the circumstellar 
environment 

In the previous analysis we have ignored the larger struc- 
ture that surrounds the circumstellar nebula. Between the 
nebula and the observer are the unshocked slow wind, the 
large bubble created by the shocked main sequence wind as 
well as a shell of interstellar matter that has been swept up 
during the main sequence phase. Beyond this shell, is the 
undisturbed ambient medium. 

The interaction between the main sequence wind and 
the ambient medium actually reduces the column density 
and helps make the circumstellar structures visible. If the 
star is embedded in a molecular cloud, the column den- 
sity between the star and the observer will make observa- 
tions impossible. However, during the main sequence phase, 
the stellar wind sweeps up the surrounding material in a 
shell, which can break out of the cloud, allowing us a direct 
view of the star and its immediate surroundings. This can 
be s een in the Rosetta Nebula around clust er NGC 2244 
(e.g. iPhebs fe Ladal 1199ft iLi fe Smith! l2QQ5h . The swept- 
up shell can contain several tens of thousands of solar 
masses. However, this actually increases the visibility, be- 
cause the swept-up shell contributes less to the column den- 
sity than if the same amount of mass had not been swept 
up and was still floating in a sphere around the star (see 



Appendix IB]) . The shocked main sequence wind itself does 
not contribute significantly to the column density. The total 
amount of mass lost during the main sequence is typically 
< 5 M and spread over a large volume. Simulations show 
main-seq uence bubbles of massive s ta rs to have radii of 
30-50 pc (iGarcia-Segura et alj ll996bllaL IFrever et all [2003, 
200 61 (Pwark adas 2005, 2007 : Ivan Marie et ■ all I2005L 1200ft 
Eldridge et al.l 120061 : iToala fc Arthurl 1201 lh . Assuming an 
ambient medium with a density of = 2/cm 3 , a main se- 
quence bubble of R = 40 pc that contains 5 M of shocked 
wind material and a swept-up shell with a cross-section of 
2pc, the column density becomes: 9.33 x 10 16 nh/cm 2 for 
the shocked main sequence wind and 8.77 x 10 19 nh/cm 2 
for the swept up shell. The column density of the un- 
shocked interstellar medium will, of course, depend on the 
distance between the observer and the star. Inside the 
shocked main sequence wind, the slow wind (RSG or LBV) 
also adds to the column density. Its contribution is diffi- 
cult to estimate, because part of it has the 1/r 2 density 
profile of the free-streaming wind and part is contained 
in a thin shell against the termination sho ck of the main 
sequence w i nd (|Garcia-Segura et al.|ll996bllal; IFrever et al.l 
20031 120061: iDwarkadasI I20Q5L 1200ft Ivan Marie et al.ll2005l 
2007; IToala fc Arthurl 1201 lh . 



The situation becomes even more complicated when 
we consider that massive stars typically occur in clusters, 
which increases the amount of matter between the neb- 
ula and the observer. Still, observations of such nebulae 
as NGC 6888 already provide as with a detailed view of 
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Fig. 17. Particle density squared [1/cm 6 ] volume rendering of the high resolution WR-LBV interaction (simulation B4) 
at the same time as the last panel of Fig. [TTJ 



the structure of the circumstellar shell. Infra-red imaging, 
which can penetrate deeper will improve on this even fur- 
ther. 



6. Discussion 

6.1. Morphology of circumstellar nebulae 

From the results shown in Section [3j it is clear that the 
morphology of a swept-up, circumstellar shell varies greatly 
depending on the wind parameters. The extremely high 
density of the slow- moving RSG wind provides far more 
resistance to the expanding shell than the faster, less 
dense, LBV wind. As a result the shell between the WR 
wind and the RSG wind is far more unstable than for 



no t ed by 

Garcia-Segura et al. (1996b, a), by Ivan Marie et all (2005, 
20071 ) and by iToala fc Arthurl ((20111). Comparison with 



these papers shows that our WR-RSG simulations show 
similar instabilities to previous work, but the RT fingers 
found in our simulations tend to be broader. This can be 
understood, when we compare the total expansion speed of 
the circumstellar shell, which for our simulations is consid- 
erably slower than what was found in earlier papers (See 
Section 16. 4ft . This lower expansion speed gives the RT fin- 
gers more time to expand sideways. Also, because we keep 
all the gas at a minimum temperature of 10 000 K, the 
thermal pressure in the hig h density areas remains high. 
iGarcia-Segura et al.l (|l996aj ) allowed the gas to cool fur- 
ther. Therefore, the high thermal pressure in the shocked 



WR wind w ould compress the RT fi ngers more than in our 
simulations. IToala fc Arthurl (|201lf ) used radiative transfer 
to determine the degree of ionization (and therefore the 
temperature) of the gas. They find that in the later stages 
of the WR-RSG interaction the temperature in the high 
density clumps tends to be at a few 1 000 K and the RT 
fingers in their models are more stretched out than ours. 
Our WR-LBV results also resemble previous simulations, 
but seem less unstable, which is probably due to the lower 
masslo ss rates in our model. For example, IToala fc Arthurl 
(2011) use models for the 60 M star that have a mass- 
loss rate reaching 10 _3 M yr _1 . The one large instability 
we find in the high-resolution 2-D model (simulation B2, 
Fig. |8j) is probably a rare event, since it is not found in 
the 3-D results. Also, they may well become less prominent 
over time as shown in Fig. [9] 



6.2. Influence of resolution 

From our simulations, it is clear that a high resolution is 
very desirable as it can significantly change the result. This 
is especially clear in the 2-D WR-LBV simulations (B1-B2), 
which show that the low resolution models do not fully re- 
solve the instabilities, which inhibits their growth. Also, the 
growth rate of the thin-shell instabilities depends on the 
thickness of the shell, as discussed in Section SJ Because 
the shell in the WR-LBV interaction moves fast, this can 
leave the instabilities with insufficient time to develop, if 
the resolution is too low to fully resolve the shell. A high 
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Iogl0(n2) Iogl0(n2) Iagl0(n2) 




Fig. 18. Frontal aspect of the circumstellar shell from Fig. fT5l (simulation A4 after 39 200 years ) showing only those parts 
with log(n 2 ) > 2 (left panel) and log(n 2 ) > 3 (centre panel) and log(n 2 ) > 4 (right panel). The larger filaments appear to 
become thinner and the smaller filaments disappear entirely. The open space between filaments seems to become much 
larger. 




Fig. 19. Frontal aspect of the circumstellar shell from Fig. [T71 (simulation B4 after 11800 years ) showing only those 
parts with log{n 2 ) > 1 (left panel), log(n 2 ) > 2 (centre panel) and log(n 2 ) > 3 (right panel). The qualitative structure 
does not change much, but the amount of "empty space" between filaments increases and the filaments themselves are 
fading. A further reduction in visibility would remove the filaments completely and leave only a handful of knots. 



resolution also allows a higher compression, which in turn 
increases the density. Since RT instabilities depend on den- 
sity contrast, this increases the growth rate. The WR-RSG 
simulation with its slow moving shell is less vulnerable, but 
the high resolution models show more structure in the in- 
stabilities than is found in low resolution. The necessity 
to use high-resolution grids is an important consideration 
when one has to decide between 2-D and 3-D simulations, 
as high resolution 3-D models will quickly become com- 
putationally expensive. Our models of the WR-LBV wind 
interaction show that the shape and size of the individual 
instabilities does not change significantly between simula- 
tion B2 and the model with twice the resolution, so the 
resolution of B2 is sufficient. 



6.3. 2-D vs. 3-D 

The main difference between 2-D and 3-D models seems to 
lie in the speed at which the instabilities develop. That this 
happens faster for 3-D models can be easily understood. If 
matter moves along the plane of the shell (as is the case for 
Vishniac instabilities) in 2-D it can only do so along the lat- 
itudinal axis, since the shell is represented as an essentially 
1-D feature. In a 3-D model, it can move along a planar 
shell structure, giving it an extra degree of freedom (the 
longitudinal axis). As a result, the local density will vary 
with the angular velocity squared, rather than as a linear 
function. As a result, at any given time, the instabilities 



will be more pronounced for a 3-D model than for its 2-D 
equivalent, as was shown by the Fourier analysis. 

The transition from 2-D to 3-D also influences the 
shape of the instabilities. For example. I Young et al.l (2001) 
showed that RT instabilities in 3-D show more structure 
than their 2-D counterparts. However, the differences are 
relatively small compared to the total instabilities and will 
likely prove impossible to observe. The behaviour of these 
instabilities is clearly different from t he convection flows 
desc ribed bylMeakin fc Arnettl ([2007D : lArnett et al.l (|2009l ) 
and lArnett fc MeakirJ (|201lh , which showed a marked dif- 
ference between the 2-D and 3-D simulations. 

It is necessary to use 3-D models if one wishes to com- 
pare the morphology of the simulated shells directly with 
observations as shown in Section El Of course, it is pos- 
sible to rotate a 2-D result around the polar axis to cre- 
ate a pseudo-3-D grid, but the result would be radically 
di fferent. This technique was use d to create projec tions 
by iGarcia-Segura et al.l (|l999f ) and IChita et al.l (|2008f ) and 
yields valuable results, but instabilities in the 2-D result 
create rings around the axis when treated in this fashion 
as can be seen in these papers. Only by actually making a 
3-D model can we demonstrate what the instabilities would 
look like. 

6.4. Expansion velocity 

It is possible to obtain the expansion speed of a circum- 
stellar nebula by measuring the blue-shift of the emis- 
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sion and/or absorption lines created by the swept-up shell. 
Because the expansion speed is a direct result of the wind 
parameters, it can be used as a quantitative means to check 
the simulation results against observational data. 

The expansion speed of the shell can be approximated 
analytically for a purely adiabatic interaction. This predicts 
an expansion velocity of 



• 2ta\ 1/3 

3M 



(6) 



with V s the bulk motion of the shell, v and m the ve- 
locity and massloss rate for the fast wind and V and M 
the ve locity and m assloss rate for the slow wind (Eq. 12.23 
from lKwok|[2QQQ[ ). This would give us an expansion speed 
of 110 kms -1 for the WR-RSG interaction and 175 kins" 1 
for the WR-LBV interaction. In our simulations we find 
85 kms -1 and 335 kms -1 respectively by measuring the 
distance travelled by the shell during the simulation. (N.B. 
In the case of the WR-RSG interaction the expansion speed 
actually varies between 80 and 89 kms -1 due to the insta- 
bilities distorting the shock. The speed of the WR-LBV 
interaction, which has a smooth forward shock, can be de- 
termined with an accuracy of about 1 kms -1 .) 

For the WR-RSG interactions, the analytical prediction 
is too high. This is partly due to the extensive instabil- 
ities in the WR-RSG simulations, which use up part of 
the available energy. Furthermore, due to the high den- 
sities involved, the WR-RSG interaction has strong radia- 
tive cooling, which allows energy to leak out of the system. 
The analytical value for the speed of the WR-LBV inter- 
action is clearly impossible, since it is actually less than 
the wind velocity of the LBV wind, so according to this 
prediction there would not be a shell. The explanation lies 
in the assump tions on wh ich the analytical model is based. 
Specifically, iKwokl ([20001 ) assumed that the mechanical en- 
ergy of the slow wind is much less than for the fast wind. 
This is correct for the AGB to post-AGB transitions that 
form planetary nebulae as well as for the WR-RSG inter- 
action. However, for the WR-LBV interaction, both winds 
have similar mechanical energy. 

The expansion velocity we find for the WR-RSG in- 
tera ction is significantly lower t h an the values found b y 
e.g. IGarcia-Segura et all (jl996aft : Ivan Marie et al.l (|2005l ): 
iFrever et al.l (|2006h for a similar situation. These mod- 
els used wind parameter s based on the massloss rates 
from lde Jager et ail ((1988), which tend to overestimate the 
massloss rates for hot stars as they ignore the effect of 
clumping in the stellar wind. Our m odels use a lower WR 
massloss rate, based on the work of iBouret et al.l (|2005l ); 
IVink fc de Koterl (|2005l ): iMokiem et al.l ([20071 ). which take 
the effect of clumping into account. Observations of WR 
wind nebulae typically give low (~ 50 — 100 kms -1 ) 
expansion veloc it ies for the WR-RSG interaction (E.g . 
Chu et all 119831: iGoudis et al.l Il988t [Smith et al.l 119881 : 



Dyson fc Ghanbari 1989; iGruendl et al. 1 12000ft . This seems 
to confirm the validity of our results. 

7. Conclusions 

In the past, 2-D hydrodynamical simulations have been 
successful in reproducing the circumstellar shells of mas- 
sive stars. We have now extended these models to 3-D. We 
find that, although the 2-D models resemble the 2-D cuts 



made through the 3-D results, the 2-D models are insuffi- 
cient to fully model the structure of circumstellar nebulae. 
Specifically, local density fluctuations in the shell that re- 
sult from hydrodynamical instabilities appear in 2-D mod- 
els as individual clumps. In reality, the instabilities form a 
lattice of filaments that can only be reproduced in a 3-D 
model. Even the slices tend to show more pronounced in- 
stabilities than the 2-D models as shown in their respective 
power spectra, though this may be nearly invisible. 

The expansion speed of the circumstellar shell in 
the WR-RSG interaction is markedly slower in ou r 
simulations than foun d by IGarcia-Segura et al.l (|l99 6a): 
Ivan Marie et al.l (|2005l ): IFrever et al.l (|2006l ). which can be 
accounted for by the lower RSG wind velocity (10 kms -1 
vs. 15 kms -1 ) and lower massloss rate during the WR 
phase (10 -5 M yr -1 vs. (10 -45 M yr -1 ). Since our 
expansion velocity (< 100 kms -1 ) lies closer to observed 
nebula expansion velocities this can serve as an additional 
validation for the lower WR massloss rates. The lower ex- 
pansion speed also has consequences for the shape of the 
instabilities, which have more time to form and expand, 
leading to broader, more irregular shapes than found in 
previous models. 

For the future we hope to expand our work to in- 
clud e stellar rotation ([Chita et al.l 120081 : Ivan Marie et al.l 
120081 ) as well as the in teraction between the winds of 
binary companions ( e.g. iPjttardl 120091 : iPittard fc Parkin! 
2010; Ivan Marie et al. 2 011b[ ) and massive stars in clusters 



([van Marie et al. 1 120121^ 7 We also intend to add additional 
physics, such as the influence of dust ([van Marie et al.l 
l2011cHCox et al.ll2012[ ). 
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Appendix A: Fourier analysis of instabilities 

We use Fourier analysis (Russl I1995L p.283-287) to quan- 
tify the instabilities in circumstellar shells (See Section [3]). 
This was done using the fast-Fourier transform (FFT) rou- 
tine provided by IDL 7.0. Calculating the Fourier trans- 
form as shown in Figs. [T2lfT3l requires the following steps. 
First of all, we map the data onto a uniform grid with 
the same resolution as the highest level of the adaptive 
mesh. This is accomplished using the same interpolation 
algorithm as used by MPI-AMRVAC for the adaptive mesh 
during the simulation. We calculate the particle density ac- 
cording to n = p/rrih- The particle density is a more prac- 
tical quantity than the mass density, because the values 
are close to one. We then isolate the part of the grid con- 
taining the circumstellar shell. In this region we calculate a 
1-D FFT along each latitudinal gridline and take the power 
spectrum, which is given by: 

P(u) = \S{v)\ 2 (A.l) 

with P(y) the power spectrum and S{v) the result of the 
FFT. We then take the average of the 1-D results and nor- 
malize it to the height of the zero-order maximum. The 
result is then plotted as shown in Figs. [T2]fT3l 



Appendix B: Column density of a swept-up shell 

As the main sequence wind of a star sweeps up the ambi- 
ent medium, the effective column density changes, because 
the matter is concentrated in a shell. As the shell expands, 
the surface over which the matter is spread out increases, 
reducing the amount of mass that a ray encounters as it 
passes through the shell. In the case of undisturbed matter 
the column density is simply p c = pism R, with p c the col- 
umn density, Pism the ambient medium density and R the 
distance to the star. If the same amount of mass has been 
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swept up in a shell with a thickness of AR that extends 
from R — AR to i?, the column density becomes 

Pc ~ ^{R-AR) 2 AR^ (BA) 
1 R 2 

^ 3 P1SM R-2AR> (B ' 2) 

for AR <C R. Therefore, the column density in the swept 
up shell is less than that of the same medium in undisturbed 
condition if 

1 R 2 

PismR > 3^ism ^_ 2AR , (B.3) 

R > 3AR. (B.4) 

For the bubbles around massive stars this is the case, so 
we can conclude that the main sequence wind reduces the 
column density by sweeping up the ambient medium. 



20 



